%% split SEDAC into countries and into connected countries in Europe
% aggregate the population data to evenly spaced points
% at the center of 0.1 degree cells
% each point is assigned the population of its cell
% population is further aggregated to 0.5 degree cells by
% make_country_graps.m

clear
close all
clc

set(0,'DefaultFigureWindowStyle','docked','DefaultFigureVisible','on')

% -----------------------
% DEFINE FOLDER LOCATIONS
folders;

%% load country list

country_list_short;

%% load Europe grid from SEDAC

% load TIF, this is a 250m cell grid
[ A,R ]  = geotiffread( [datafolder_SEDAC,....
                        'gpw-v4-population-count-adjusted-to-2015-unwpp-country-totals_2005.tif'] );
                   
% zeros appear as tiny negative numbers but they add up, set them to zero                    
A( A<0 )=0;

% Locate NW corners of each cell

% Y runs North-South
Ygrid = R.LatitudeLimits(2):-R.CellExtentInLatitude:R.LatitudeLimits(1)+R.CellExtentInLatitude;

% X runs West-East
Xgrid = R.LongitudeLimits(1):R.CellExtentInLongitude:R.LongitudeLimits(2)-R.CellExtentInLongitude;

%% aggregate into coarser cells and generate places map structure

% set cell size, in minutes
cellsize = 0.1;
NN = round( cellsize/R.CellExtentInLatitude ); % each cell in coarser grid will contain NN^2 cells in finer original grid

% centers of cells
Ygrid_coarse = max(Ygrid)-cellsize/2:-cellsize:min(Ygrid);
Xgrid_coarse = min(Xgrid)+cellsize/2:cellsize:max(Xgrid);

clc
counter = 1;
for j=1:length(Xgrid_coarse)
    
    [ j length(Xgrid_coarse) ]
    
    for i=1:length(Ygrid_coarse)
        
        places(counter).Geometry = 'Point';
        places(counter).X = Xgrid_coarse(j);
        places(counter).Y = Ygrid_coarse(i);
        
        temp = A( NN*( i-1 )+1:NN*i,...
                  NN*( j-1 )+1:NN*j );
              
        places(counter).population = sum( temp(:) );  % the center of each cell is assigned its population
        
        counter = counter+1;
        
    end
end
                    
%% retain the section of grid corresponding to each country

for country_n=1:Ncountries
    
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name

    if ~ispc()   %% we are in a mac

        % folder to save the outcome
        datafolder_LAC_country = [ datafolder_LAC,country_icc,'/'];

    else   %% we are in a pc   

        % folder to save the outcome
        datafolder_LAC_country = [ datafolder_LAC,country_icc,'\' ];

    end
    
    clc
    disp( ['Country: ',country] )    
    
    % load country_bounds
    file  = [ datafolder_LAC_country,'country_bounds.mat']; 
    load(file);
    
    % keep places within the defined boundaries
    keep_pop = inpolygon( cell2mat({places.X}'), cell2mat({places.Y}'), country_bounds.X, country_bounds.Y );
    places_sedac = places( keep_pop );
    
    % save place
    file  = [ datafolder_LAC_country,'popplaces_sedac.mat'];
    save(file,'places_sedac');
    
    %% plot
    
    h=figure;
    mapshow(places_sedac,'Marker','.','MarkerSize',10);
    mapshow(country_bounds,'FaceColor','none')
    totpop =sum( cell2mat( {places_sedac.population} ) ) ;
    title([country,'. Tot Pop (Million): ',num2str(round(totpop/10000)/100)]); 
    if ismac
        saveas(h, [datafolder_SEDAC,'maps/',country,'_map_sedac.pdf']) 
    else
        saveas(h, [datafolder_SEDAC,'maps\',country,'_map_sedac.pdf']) 
    end
        
end

